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INTRODUCTION 


Millions  of  Americans,  including  a  large  percentage  of  active  military  personnel  and  veterans,  are  affected  by 
three  major  blinding  diseases  of  the  retina  and  optic  nerve:  diabetic  retinopathy,  age-related  macular 
degeneration,  and  glaucoma.  These  diseases  typically  affect  the  elderly,  but  a  significant  number  of  young 
people  also  fall  victim.  As  a  result  of  current  advances  in  surgical  and  pharmacologic  therapeutics,  vision  loss 
due  to  all  these  diseases  can  be  halted  or  even  reversed  if  the  disease  is  detected  early.  New  advanced  detection 
methods  are  available,  but  are  only  interpretable  by  very  experienced  specialists.  The  goal  of  this  research  is  to 
advance  the  frontiers  of  a  promising  imaging  technology  called  optical  coherence  tomography  (OCT), 
developing  advanced  instrumentation  and  new  software  applications  to  enable  high  resolution  evaluations  of  the 
retina  and  optic  nerve  that  could  be  easily  performed  by  non-physician  health  care  professionals.  In  addition  to 
on-site  use,  such  techniques  could  be  further  applied  for  telemedicine  transmittal  of  images  (from  the  battlefield 
or  other  remote  sites)  for  further  evaluation  by  specialist  ophthalmologists.  The  hypothesis  to  be  tested  is  that 
advancements  in  OCT  technology  will  enhance  military  and  nonmilitary  ocular  health  capabilities.  The  three 
specific  research  objectives  are  as  follows:  1)  Develop  methods  for  accurate  vertical  alignment  and  matching, 
including  the  overlay  of  a  color  or  red-free,  fundus  photograph  on  an  OCT  image  for  registration.  2)  Develop 
the  mathematical  algorithms  to  analyze  and  quantify  OCT  datasets  automatically.  3)  Further  improve  the 
resolution  of  the  spectral  OCT,  correlate  images  with  pathology,  and  develop  algorithms  for  3-D  visualization 
of  datasets. 

BODY 

Task  1.  Development  of  algorithms  for  quantitative  evaluation  of  retinal  pathology  and  for  monitoring 
disease  progression  and  response  to  treatment. 

a.  Develop  3D  segmentation  algorithm  for  quantitative  mapping  of  the  retinal  thickness  using  current  3D  OCT 
data  for  normal  and  diseased  human  eyes. 

The  development  of  high-speed  spectral-domain  OCT  systems  allows  for  the  introduction  of  new  imaging 
modalities,  including  the  acquisition  of  three-dimensional  datasets.  Detailed  information  about  the  retinal 
structure  over  large  areas  is  encoded  in  these  3-D  datasets.  A  crucial  challenge  to  exploiting  the  full  potential  of 
retinal  OCT  imaging  is  the  ability  of  extracting  reliable,  quantitative  information  from  the  scans.  Our  efforts 
towards  providing  an  answer  to  this  question  are  an  important  component  of  the  research  supported  by  the 
grant.  In  particular  we  developed,  implemented  and  validated  3D  segmentation  algorithms  to  map  the  retinal 
geometry  from  the  OCT  datasets.  One  of  the  main  products  of  the  segmentation  is  the  measurement  of  retinal 
thickness.  Such  measurements  provide  the  framework  for  the  quantitative  evaluation  of  retinal  pathology. 
Multimodal  registration  between  OCT  fundus  reconstructions  and  other  en  face  retinal  images  (i.e.  fundus 
photography,  fluorescein  angiography,  fundus  autofluorescence)  gives  us  the  tools  to  map  the  precise  retinal  site 
corresponding  to  each  image  pixel  and  track  a  specific  location  over  time  (across  successive  datasets).  This 
provides  an  unprecedented  ability  for  monitoring  disease  progression  and  response  to  treatment.  This  type  of 
information  could  prove  very  valuable  to  the  physician  in  formulating  clinical  decision  as  well  as  for  improving 
the  understanding  of  pathological  processes. 

The  first  step  involved  developing  image  processing  software  to  perform  tasks  such  as  feature  recognition  and 
edge  segmentation.  There  are  several  issues  that  make  segmentation  of  OCT  images  very  difficult  in  general. 
Some  of  the  main  complications  that  one  needs  to  overcome  are  the  presence  of  speckle  noise,  the  relative  low 
contrast  and  signal  to  noise  ratio  in  a  OCT  image  compared  to  x-ray  CT  and  MRI  imaging,  and  the  large 
variability  of  retinal  features’  appearance,  especially  in  the  presence  of  pathology.  The  difficulty  in  extracting 
quantitative  information  from  OCT  images  is  well  illustrated  by  the  problems  with  the  commercial  OCT 
instrument  (StratusOCT,  Carl  Zeiss),  whose  retinal  thickness  analysis  has  been  shown  to  generate  a  large 
number  of  artifacts.  An  additional  problem  in  3D  spectral-domain  OCT  imaging  is  the  very  big  size  of  the 
datasets.  The  very  large  number  of  A-scans  to  be  analyzed  (>40000)  greatly  restricts  the  acceptable  rate  of 
failure  of  the  algorithms  to  be  used,  before  errors  become  clearly  evident.  Also  the  computational  power/time 
needs  of  potential  algorithms  need  to  be  taken  into  careful  consideration. 
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The  algorithms  developed  by  our  group  allow  us  to  reconstruct  the  three  dimensional  geometry  of  the  patient’s 
retina.  In  particular  we  can  visualize  the  geometry  of  the  internal  limiting  membrane  (ILM)  and  the  retinal 
pigment  epithelium  (RPE)  as  surfaces  in  a  three  dimensional  space.  In  Fig.  1,  we  show  this  construct  for  the 
retina  of  a  patient  suffering  from  non-exudative  age-related  macular  degeneration  (AMD).  Retinal  thickness 
maps,  i.e.  the  distances  between  ILM  and  RPE,  can  be  then  computed  and  visualized.  The  retinal  thickness 
maps  we  can  generate  from  SD-OCT  images  are  vastly  more  detailed,  accurate,  and  reproducible  than  those 
available  with  StratusOCT  or  any  other  current  imaging  technology. 


Spectral-domain  OC Is  can  generate  3D  datasets 
covering  large  retinal  areas 


r  .  A  typical  scan  consists  of  a  raster  of  equally  spaced 

OCT  fundus  reconstruction.  It  allows  scan  B-scans  covering  a  square  region, 

registration  as  well  as  quality  assessment. 


Segmentation  algorithms  can  be  used  to 
reconstruct  the  retinal  geometry  and 
create  thickness  maps. 


ILM 


Total  Retinal  Thickness 


RPE 


Retinal  Thickness  Map 


Fig.  1  Our  analysis  of  OCT  datasets. 


The  algorithms  we  developed  are  based  on  an  iterative  architecture,  where  an  initial  guess  is  successively 
evaluated  and  improved  upon,  according  to  a  set  of  principles  that  may  be  tailored  to  the  specific  structure  to  be 
analyzed  and/or  the  particulars  of  the  disease  model  of  interest.  The  general  design  of  the  algorithms  is 
described  in  the  following  flowchart  (Fig.  2). 

We  will  assume  here  that  we  are  dealing  with  an  OCT  dataset  acquired  using  a  raster  scan  consisting  of  dA  x  dH 
x  dv  data  points,  where  dA  is  the  A-scan  size,  dH  is  the  number  of  scans  in  each  B-scans,  and  dy  is  the  number  of 
B-scans  in  the  raster  scan.  The  OCT  dataset  is  then  a  dA*  dHx  dv  dimensional  array  A(i,j,k),  where  the  pixel 
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indexes  are  equivalent  to  a  system  of  3-dimensional  Euclidean  orthogonal  coordinates.  The  main  output  consists 
of  two  dH><dv  arrays  representing  the  ILM  and  the  RPE  positions,  i.e.  the  ILM(j,k)  entry  is  an  integer  giving  the 
pixel  number  in  the  A-scan  A(.,j,k)  corresponding  to  the  ILM  leading  edge. 


Fig.  2.  Basic  Structure  of  the  Algorithm. 

The  first  step  consists  of  some  preprocessing  the  OCT  data,  whose  main  purpose  is  to  reduce  the  amount  of 
speckle  noise  and  then  computing  some  basic  statistical  characteristics  of  the  remaining  noise  process.  Also 
some  signal-level  normalization  is  achieved. 

At  this  point  a  quick  analysis  of  each  A-scan  is  carried  out  to  produce  an  initial  estimate  for  the  ILM  and  RPE 
arrays.  While  it  would  of  course  be  desirable  to  achieve  the  correct  segmentation  already  at  this  stage,  our 
approach  only  requires  that  this  initial  estimate  be  “mostly  correct”.  This  is  actually  the  main  assumption  behind 
the  design  of  the  algorithm:  if  at  least  a  sizable  percentage  on  the  boundaries  can  be  quickly  estimated  in  a 
correct  manner  (and  if  we  can  recognize  this  correct  region),  then  we  can  focus  our  computational  resources 
effectively  on  the  “bad  regions”  to  improve  the  result. 

An  iteration  loop  is  therefore  set  up  where  the  algorithm  repeatedly  improves  the  segmentation  results  by 
selecting  “bad  regions”  where  the  current  segmentation  is  “suspect”  and  analyzing  these  bad  regions  more 
carefully.  The  process  stops  when  either  there  are  no  bad  regions  left  or  the  algorithm  has  traversed  a  preset 
maximal  number  of  loops.  The  routines  which  decompose  a  given  segmentation  in  good  and  bad  regions,  and 
the  routines  that  analyze  and  correct  the  segmentation  on  the  bad  regions,  are  the  computational  core  of  the 
algorithm.  They  strive  to  encode  a  series  of  hypotheses  as  well  as  a  fair  amount  of  empirical  knowledge  on  the 
anatomy  and  geometry  of  the  retina.  The  concentration  of  this  knowledge  into  a  suitable  set  of  mathematical 
conditions  which  do  a  good  job  of  discriminating  problematic  areas  in  the  segmentation,  and  work  well  across  a 


6 


very  wide  spectrum  of  retinal  conditions,  is  the  crucial  task  to  be  achieved.  We  believe  our  choices  can  be 
shown  to  perform  quite  successfully. 


Fig.  3.  Segmentation  result  for  a  human  retina  with  RPE  detachment.  Left:  3D  view  of  the 
segmented  ILM  and  RPE;  Center:  3D  view  of  the  segmented  RPE;  Right:  the  calculated  retinal 
thickness  map. 


We  have  focused  so  far  on  tools  that  describe  the  ILM  and  the  RPE,  i.e.  the  inner  and  outer  boundaries  of  the 
retina.  However,  these  algorithms  could  be  generalized  to  provide  quantitative  information  about  the  damage  to 
different  biological  layers  caused  by  the  diseases,  like  for  instance  the  nerve  fiber  layer  (NFL),  the  ganglion  cell 
layer  (GCL),  or  the  photoreceptors.  The  ability  to  accurately  measure  such  retinal  structures  would  likely 
further  our  understanding  of  the  progression  of  various  retinal  diseases. 


b.  Develop  3D  registration  algorithm  for  monitoring  disease  progression  and  response  to  treatment  using 
current  3D  OCT  data. 


Multimodal  and  intra-modal  retinal  image  registration  is  important  to  clinical  diagnosis  and  treatment. 
Currently,  there  are  several  types  of  retinal  imaging  modalities,  including  OCT,  fluorescein  angiogram,  fundus 
autofluorescence  imaging,  and  color  fundus  photography.  Different  types  of  imaging  reveal  different 
characteristics  the  retina  by  using  different  contrast  mechanisms,  which  provide  complimentary  localized 
information  of  the  retina.  At  the  same  time,  retinal  images  taken  at  different  time  can  reveal  the  progression  and 
the  response  to  treatment  of  retinal  diseases.  We  have  developed  an  effective  algorithm  to  register  the  various 
types  of  fundus  images:  OCT  fundus  images  taken  at  different  times  and  OCT  fundus  images  to  color  fundus 
photograph  of  the  same  eye. 

Registration  methods  usually  consist  of  several  steps  [1-3]:  feature  detection;  transform  model  estimation; 
optimization  function  design;  and  optimization  strategies.  We  do  not  choose  intensity-based  feature  for 
multimodal  image  registration  because  even  for  intra-modal  retinal  images,  intensity  can  change  due  to  changes 
in  imaging  parameters.  For  retinal  image  registration,  blood  vessel  patterns  [4,5]  can  be  taken  as  a  relatively 
stable  feature.  Some  papers  use  the  vascular  landmarks  [6,7],  like  bifurcations  and  crossovers  [8].  However,  for 
low  quality  images  or  images  with  disease,  it  is  not  easy  to  detect  bifurcations  and  crossovers. 

Here,  we  use  blood  vessel  ridges  as  a  specific  feature  [9].  Defined  as  points  where  the  image  has  an  extremum 
in  the  direction  of  the  largest  surface  curvature  [10],  ridges  are  a  natural  feature  of  blood  vessels  and  are  usually 
approximately  center  lines  of  blood  vessels.  By  using  ridges  we  can  avoid  the  problem  of  determining  the  whole 
blood  vessel  area,  which  not  only  saves  computation  time  but  also  avoids  the  potential  bad  effect  of  imprecisely 
detected  vessel  area  on  registration. 

Among  different  transformation  models,  like  translation,  affine,  and  quadratic  models  [5,  6],  affine 
transformation  is  relatively  simple  while  generating  acceptable  results: 
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where,  the  point  ( x,y )  in  one  image  will  be  transformed  to  (x',y')  in  the  other  image  by  the  affine 


transformation  matrix 
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,  which  includes  the  combination  of  translation,  rotation,  shearing  and 


scaling.  In  the  transformation  matrix,  a3l  and  ai2  are  translation  parameters,  an,  al2,  a2l,  and  a22  are 

parameters  related  to  shear,  rotation  and  scale.  We  use  the  affine  transformation  model  in  our  retinal  image 
registration  algorithm. 


Optimization  function  is  to  evaluate  the  similarity  between  a  pair  of  images.  For  intensity  based  matching  [1], 
correlation,  mutual  information,  and  Fourier  transformation  can  be  used  as  an  optimization  function.  For 
registration  with  points  as  features  [7,  5],  distance  is  a  good  similarity  measurement.  In  literature  there  are  a  lot 
of  similarity  functions  [11-13]  to  compute  distance  between  points. 


After  the  optimization  function  is  chosen,  we  need  to  choose  a  proper  optimization  strategy  [2],  Here,  we  focus 
on  optimization  strategies  using  points  as  features,  since  our  algorithm  uses  vessel  ridges  as  features.  Ref.  6 
used  a  hierarchy  of  models,  a  random  sampling  search  technique  and  iteratively-reweighted  least-squares  to 
estimate  the  quadratic  transformation.  A  dual-bootstrap  iterative  closest  point  algorithm  is  used  in  Ref.  7. 
Usually  initialization  does  not  yield  points  suitable  for  global  optimization  except  for  some  local  area,  so  it  is 
not  necessary  to  match  all  the  points.  Since  the  alignment  in  the  small  initial  region  is  reasonably  accurate,  it 
helps  to  expand  registration  area  gradually.  Ref.  5  used  brute  force  search  to  get  translation  estimation,  then 
estimated  the  affine  or  quadratic  transformations  by  iterative  closes  point  algorithm.  The  test  fundus  images  are 
at  a  similar  scale. 


Our  algorithm  adopts  vessel  ridges  as  features.  For  different  fundus  images,  they  have  different  resolutions. 
According  to  the  resolution  parameters,  we  rescale  fundus  vessel  ridge  images  to  a  similar  scale.  Brute  force 
search  is  first  used  to  estimate  the  translation  parameters  and  scale  parameter.  Then  the  translation  and  scale 
estimation  with  the  largest  similarity  is  taken  as  the  initialization  of  the  iterative  closes  point  algorithm  (ICP)  to 
get  a  more  accurate  transform.  Experiments  show  that  affine  transform  works  well. 

Three  steps  were  taken  in  the  registration  algorithm: 

Step  1:  Vessel  ridges  are  detected  for  different  fundus  images. 

Step  2:  Brute  force  search  is  used  to  estimate  the  translation  parameters  and  scale  parameter.  For  different  types 
of  fundus  images,  rescale  parameter  and  the  size  of  the  matching  area  are  estimated  as  input  parameters. 

Step  3:  With  the  estimation  of  translation  and  scale  parameters,  the  iterative  closest  point  algorithm  (ICP)  is 
used  to  find  a  more  exact  affine  transform. 


Fig.  4  test  registration  result  for  OCT  fundus  image  vs  OCT  fundus  image.  The  images  were 
taken  separately  in  the  same  day  for  a  patient,  a)  base  image,  (b)  test  image,  (c)  the 
registration  of  detected  blood  vessel  ridges  for  the  two  images,  (d)  registered  test  image. 

shows  the  registration  result  for  an  OCT  fundus  image  vs.  OCT  fundus  image.  The  OCT  images  were  acquired 
with  Carl  Zeiss’  Cirus  separately  in  the  same  day.  Fig.  5  shows  the  test  result  for  another  patient  imaged  at 
different  dates.  Fig.  6  shows  the  comparison  of  registration  algorithms  based  on  intensity  with  that  based  on 
blood  vessel  ridges.  We  can  see  that  the  feature  of  vessel  ridges  is  more  stable  than  that  of  intensity.  We  used 
Image J  software  in  Ref.  14  to  register  the  two  examples  for  comparison  because  our  registration  algorithm  can 
not  be  applied  to  intensity  based  registration.  Fig.  7  shows  the  test  result  fro  registration  of  the  OCT  fundus 
image  on  the  color  fundus  photo. 
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(c)  The  registration  of  detected  blood  (A)  Registration  result, 

vessel  ridges  for  the  two  images. 

Fig.  4  test  registration  result  for  OCT  fundus  image  vs  OCT  fundus  image.  The  images  were 
taken  separately  in  the  same  day  for  a  patient,  a)  base  image,  (b)  test  image,  (c)  the 
registration  of  detected  blood  vessel  ridges  for  the  two  images,  (d)  registered  test  image. 
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(a)  Base  image 


(b)  Test  image 


(c)  The  registration  of  detected  blood  (d)  Registration  result, 

vessel  ridges  for  the  two  images. 


Fig.  5  test  registration  result  for  OCT  fundus  image  vs  OCT  fundus  image.  The  images  were 
taken  at  different  dates. 
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(a)  Base  image 


(b)  The  base  image  is  under  the  registered 
test  image  based  on  intensity  by  ImageJ 
software.  We  can  see  clearly  the  deviation 
of  corresponding  vessels. 


(c)  The  base  image  is  under  the  registered 
test  image  based  on  vessel  ridges  by 
ImageJ  software. 

Fig.  6  Comparison  of  registration  algorithms  based  on  intensity  with  that  based  on  blood 
vessel  ridges.  We  can  see  that  the  feature  of  vessel  ridges  is  more  stable  than  that  of  intensity. 
We  used  ImageJ  software  in  Ref.  14  to  register  the  two  examples  for  comparison  because  our 
registration  algorithm  can  not  be  applied  to  intensity  based  registration. 
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(a)  The  registration  result.  The  OCT  fundus  image  is 
superimposed  on  the  color  fundus  image. 


(b)  The  registration  of  detected  blood 
vessel  ridges  for  the  two  images. 


Fig.  7  test  registration  result  for  OCT  fundus  image  vs  color  fundus  photo. 


c.  Test  the  accuracy  of  the  segmentation  algorithm  by  comparing  the  segmentation  results  done  by  the 
algorithm  and  by  experts,  respectively. 

The  performance  of  a  segmentation  algorithm  like  the  one  described  here  can  be  judged  in  terms  of  two 
important  measures:  accuracy  and  reproducibility.  These  characteristics  are  not  necessarily  related  and  are  both 
important  for  clinical  applications.  Also  algorithms’  performance  should  be  evaluated  on  the  full  range  of 
pathologies  presented  in  the  patient  population.  The  accuracy  of  any  new  algorithm  has  to  be  tested  by 
comparing  the  retinal  boundaries  generated  by  the  algorithm  with  boundaries  manually  generated  by  an  expert, 
e.g.  a  clinician  or  a  pathologist  who  is  familiar  with  OCT  image  interpretation.  This  information  can  be  used  to 
evaluate  the  performance  of  the  algorithms,  as  well  as  to  quantify  further  improvements.  Testing  the 
reproducibility  involves  scanning  a  given  set  of  eyes  repeatedly.  Because  it  is  often  difficult  to  image  exactly 
the  same  retinal  area  every  time,  the  ability  to  perform  image  registration  becomes  an  important  component  of 
assessing  the  reproducibility  of  measurements  in  the  context  of  spectral  domain  OCT. 

For  the  analysis  of  the  segmentation  algorithm’s  performance  we  focused  on  raster  scans  covering  a  6x6x2  mm 
volume  and  acquiring  either  200x200  A-scans  equally  spaced  on  the  retina  or  512x128  A-scans  covering  the 
same  region.  Automated  segmentation  of  the  internal  limiting  membrane  (ILM)  and  anterior  retinal  pigment 
epithelium  (RPE)  boundaries  generates  surfaces  in  3-D  space  and  retinal  thickness  maps. 

132  eyes  from  patients  imaged  at  Bascom  Palmer  were  randomly  selected  representing  the  full  range  of  retinal 
diseases.  Also  included  were  12  eyes  of  normal  volunteers.  For  each  eye  a  SD-OCT  dataset  was  entered  in  the 
study.  The  OCT  datasets  were  divided  amongst  three  retina  specialist,  who  traced  manually  the  ILM  and  RPE 
boundaries  a  number  of  predetermined  B-scans.  These  manually  drawn  boundaries  were  then  compared 
pointwise  with  the  results  of  the  automated  segmentation. 

The  statistical  distributions  of  the  differences  (pixelwise)  between  the  manual  and  automated  segmentations  are 
described  separately  for  ILM  and  RPE,  for  both  normal  eyes  and  eyes  with  disease.  In  addition  to  these  error 
probability  distributions,  we  defined  a  threshold  for  registration  failure.  A  computer  generated  boundary  was 
considered  a  failure  if  its  distance  from  the  corresponding  “true  edge”  was  above  20  pixels  at  least  10%  of  the 
time.  Of  course  any  such  failure  criteria  are  somewhat  arbitrary,  nevertheless  we  believe  it  come  close  to 


12 


capturing  a  clinical  judgment  about  the  usefulness/appropriateness  of  a  given  segmentation  result.  Two 
particular  examples  of  algorithm  failure  are  shown  in  Fig.  8. 


Fig.  8  Two  B-scans  showing  failures  of  the  RPE  segmentation.  Failures  can  be  fairly  subtle  (example  on 
the  left)  in  the  example  on  the  right,  or  sometime  substantial. 

Normal  subjects  are  of  course  the  simplest  case.  The  anatomical  concepts  are  very  clear,  image  quality  tends  to 
be  good,  and  segmentation  in  this  case  is  relatively  easy.  The  segmentation  algorithm  did  not  fail  on  a  set  of  36 
B-scans  from  12  eyes.  The  difference  probability  distributions,  shown  below  in  Fig.  9,  are  very  tight  with  means 
close  to  zero.  The  standard  deviation  is  somewhat  higher  for  the  RPE  than  for  the  ILM,  reflecting  the  fact  that 
identifying  the  ILM  boundary  is  in  general  a  more  straightforward  task. 


probability  distributions 
of  the  difference 
between  manual  and 
automated 
segmentation 


x-axis:  distance  in  pixels. 

A  positive  distance  correspond  to  the 
manual  boundaries  above  the  computer 
generated  one. 


RPE:  n=36  mean=  1.034.  sd=4.0SD 


Fig.  9  Accuracy  results  for  normal  eyes. 

Our  sample  of  patients’  eyes  is  meant  to  be  representative  of  what  is  typically  imaged  at  the  Bascom  Palmer 
retinal  clinic.  These  are  randomly  chosen  eyes  presenting  a  full  range  of  pathologies.  No  special  provision  was 
made  to  eliminate  scans  with  less  than  optimal  image  quality.  The  idea  was  to  obtain  results  which  would 
evaluate  the  performance  of  the  algorithm  in  a  setting  similar  to  what  clinicians  would  actually  see  in  a  busy 
tertiary  care  center.  While  the  presence  of  pathology  clearly  affects  the  performance  of  automated  segmentation, 
the  probability  distributions  remained  quite  narrow  in  all  cases. 
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ILM:  n=l45  (512x120  scans) 


ILM:  n-305 


72  200x200  datasets. 
In  red  is  the  best  fit  of  a 
gaussian  distribution. 


ILM  :  8  failures  out  of  160 
RIPE:  32  failures  out  of  212 


60  512x128  datasets. 
In  red  is  the  best  fit  of  a 
gaussian  distribution. 


ILM:  8  failures  out  of  145 
RPE:  27  failures  out  of  180 


All  (132)  patient’s  datasets 
(do  not  include  normals). 
In  red  is  the  best  fit  of  a 
gaussian  distribution. 


ILM:  mean=-1 .42,  sd=6.15 
RPE:  mean=2.56,  sd=12.64 


RPE:  m-212  (200x200  scans) 


Fig.  10  Accuracy  results  for  eyes  with  disease. 

The  variance  of  distribution  of  differences  is  larger  by  a  factor  roughly  equal  to  three.  The  segmentation  of  the 
RPE  shows  a  tendency  to  be  slightly  biased  towards  positive  errors  (i.e.  the  automated  segmentation  is  more 
often  below  the  manually  drawn  curve).  On  both  scan  patterns  the  ILM  segmentation  failed  on  ~  5%  of  the  13- 
scans,  while  the  RPE  failed  on  ~  15%  of  the  B-scans.  It  is  important  to  keep  in  mind  that  it  is  often  difficult  to 
identify  retinal  structures  (in  particular  the  “RPE”)  in  the  presence  of  pathologies  (as  a  matter  of  fact  it  is  often 
misleading  to  even  use  the  word  RPE  to  describe  what  is  really  an  outer  retinal  edge.  In  many  cases  the  actual 
RPE  can  be  compromised  and/or  missing).  Even  the  retinal  experts  could  not  sometime  agree  to  a  “right” 
position  for  the  outer  retinal  edge.  Therefore  the  results  shown  in  Fig.  10  are  quite  encouraging. 

The  reproducibility  of  the  retinal  thickness  measurements  is  estimated  by  computing  the  standard  deviation, 
both  pointwise  and  averaged  over  regions,  of  registered  thickness  maps.  Variations  in  OCT  measurements  are  in 
general  due  to  several  factors.  Datasets  acquired  with  the  new  SD-OCT  instruments  can  in  principle  be 
registered  using  the  OCT  fundus  images,  therefore  minimizing  the  component  of  variation  due  to  eye 
movements  as  well  as  scan  aiming/patient  fixation  problems. 
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Fig.  11  Registration  of  OCT  images  can  be  needed  to  analyze  reproducibility:  the  left  eye  of  a  patient 
was  imaged  three  times  over  two  months.  Clearly  the  scans  do  not  overlap  exactly. 

Reproducibility  is  excellent  in  normal  subjects.  In  this  case  the  retinal  thickness  changes  continuously  and  the 
thickness  maps  can  be  registered  in  an  effective  manner.  We  see  in  Fig.  12  a  typical  example  where  five 
separate  OCT  datasets  of  a  normal  eye  were  acquired  and  compared  after  registration.  The  pointwise  standard 
deviation  map  shows  that  at  most  pixels  the  standard  deviation  of  the  five  thickness  measurements  is  well  below 
5pm  (and  always  below  8pm).  The  average  standard  deviation  is  2.52  pm. 
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5  separate  scans  of  the  same  normal  eye 


Pointwise  standard  deviation  map 
Average  pointwise  sd=2.52  pm 


Fig.  12  Reproducibility  for  a  normal  eye. 
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Pointwise  standard  deviation  map 
Average  pointwise  sd=6.67  pm 


Fig.  13  Reproducibility  for  an  eye  with  AMD. 


Fig.  9. 
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In  the  presence  of  pathologies  the  reproducibility  of  thickness  measurements  might  deteriorate  somewhat.  The 
sharp,  localized  changes  in  thickness  maps  due  to  drusen,  for  instance,  can  increase  the  standard  deviation  as  we 
see  in  Fig.  13.  This  effect  may  be  partly  due  to  difficulties  with  the  registration,  which  in  this  case  is  particularly 
important  and  may  not  be  as  good  as  one  would  like.  It  should  be  noted  that  both  the  average  thickness  map  and 
the  standard  deviation  map  show  some  black  area  along  the  left  and  lower  edges.  These  are  the  areas  where  the 
three  dataset  fail  to  overlap.  It  might  also  be  remarked  that,  while  the  pointwise  standard  deviation  map  show 
small  areas  with  relatively  large  standard  deviation  values,  the  qualitative  look  of  the  three  thickness  maps  is 
very  close.  For  most  clinical  purposes  average  retinal  thickness  values  over  regions  of  interest  might  be  the 
relevant  measurement.  For  instance  physicians  are  familiar  with  the  ETDRS  grid  which  defines  nine  regions 
bounded  by  three  concentric  circles  (Fig.  10).  When  using  the  ETDRS  grid  to  produce  average  retinal  thickness 
values,  the  scans  in  Fig.  9  show  a  very  high  degree  of  reproducibility. 


Average  retinal  thickness  Standard  deviation  in 
in  each  sector  each  sector 

Fig.  14  Average  thickness  and  standard  deviation  over  the  9  ETDRS  regions  for  the  eye  in  Fig.  13. 

The  3D  datasets  generated  by  our  SD-OCT  instrument  can  be  automatically  analyzed  with  the  segmentation 
software  we  developed.  This  allows  us  to  generate  thickness  maps  consisting  of  40000-70000  data  points 
distributed  over  a  square  region  of  the  retina.  This  wealth  of  sample  points  creates  images  that  can  accurately 
describe  even  small  features. 

The  boundaries  selected  from  the  automated  algorithm  are  generally  in  very  good  agreement  point-wise  with 
the  manually  drawn  boundaries.  In  particular  on  diseased  eyes  these  lines  are  within  20pm  of  each  other  on  well 
over  90%  of  pixels  about  85%  of  the  times.  It  should  be  kept  in  mind  that  many  of  these  eyes  were  seriously 
diseased,  poor  image  quality  was  not  an  exclusion  criteria,  and  different  retinal  specialists  were  often  at  odds  on 
where  boundaries  (especially  the  RPE)  should  be  drawn.  The  standard  deviation  of  the  probability  distribution 
for  the  errors  was  below  10  pm,  when  controlling  for  outliers. 

The  variance  of  the  segmentation  of  different  scans  of  the  same  eye  was  measured  at  every  pixel  after  the 
images  were  registered.  This  point-wise  variance,  as  well  as  the  mean  variance,  was  typically  well  below  5  pm. 

Task  2.  Develop  ultra-high  resolution  ophthalmic  OCT  system. 

A  high-speed  high  resolution  3D  SD-OCT  was  built.  A  schematic  of  the  experimental  system  for  the 
preliminary  studies  is  shown  in  Fig.  15.  In  the  SD-OCT  system,  the  low-coherence  light  from  a  three-module 
superluminescent  diode  (T-840  Broadlighter,  Superlum  Diodes  Ltd.  Moscow,  Russia)  with  center  wavelength 
of  840  nm  and  FWHM  bandwidth  of  100  nm  is  coupled  into  an  optical  fiber-based  Michelson  interferometer.  In 
the  sample  arm,  the  sample  light  is  delivered  to  the  retina  by  a  modified  optical  head  of  an  OCT  2  system  (Carl 
Zeiss  Meditec  Inc.,  Dublin,  CA).  The  power  of  the  sample  light  was  lowered  to  750pW  by  adjusting  the  source 
power  to  ensure  that  the  light  intensity  delivered  to  the  eye  was  within  the  ANSI  standard.  In  the  detection  arm, 
a  spectrometer  consisting  of  a  collimating  lens,  a  transmission  grating  (1200  line/mm),  a  multi-element  imaging 
lens  (f  =  180  mm),  and  a  line  scan  CCD  camera  (Aviiva-SM2-CL-2014,  2048  pixels  with  14  micron  pixel  size 
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operating  in  12-bit  mode)  was  used  to  detect  the  combined  reference  and  sample  light.  The  calculated  spectral 
resolution  of  the  spectrometer  was  0.055  nm,  which  corresponds  to  a  detectable  imaging  depth  range  of  3.1  mm 
in  air.  OCT  scans  consisting  of  a  total  number  of  65536  depth  scans  (A-scans)  takes  2.7  seconds.  An  image 
acquisition  board  (NI  IMAQ  PCI  1428)  acquired  the  image  captured  by  the  camera  and  transferred  it  to  a 
computer  workstation  (IBM  Intelli Station  Z  Pro,  dual  3.6  GHz  processor,  3  GB  memory)  for  signal  processing 
and  image  display.  A  complete  raster  scan  consisting  of  128  x  512  scanning  steps  took  about  2.7  seconds  when 
the  A-line  rate  of  the  OCT  system  was  set  to  be  24  kHz.  At  this  operating  condition,  the  measured  sensitivity 
was  about  95dB. 


Bk7  glass 

Isolator  block 


Fig.  15.  Schematic  of  the  preliminary  experimental  system.  SLD:  superluminescent  diode 
(Broadlighter,  Superlum  Diodes  Ltd,  Moscow,  Russia);  PC:  polarization  controller.  The 
spectrometer  consists  of  a  f=50  mm  collimating  lens,  a  1200  line/mm  transmission  grating  and  a 
f=200  mm  imaging  lens.  The  CCD  camera  is  a  2048  element  linear  array  with  14  pm  element  size 
(Aviiva-M2-CL-2014,  Atmel,  USA). 

Calculation  of  the  OCT  signal 

In  spectral-domain  OCT,  the  combined  back-reflected  sample  and  reference  light  in  a  Michelson  interferometer 
is  detected  by  a  spectrometer  together  with  an  array  detector  (usually  a  CCD  camera).  The  following  analyses 
assume  that  a  linear  array  detector  is  used  whose  elements  are  aligned  in  the  direction  along  which  the  spectrum 
is  spread  in  the  spectrometer.  We  also  ignore  the  polarization  effects  in  the  interference  between  the  reference 
and  sample  light  without  losing  the  generality  of  the  analysis.  The  signal  detected  by  the  array  detector  is  called 
a  spectral-domain  signal  to  distinguish  it  from  the  time-varying  signal  detected  in  a  conventional  time-domain 
OCT.  The  light  intensity  incident  on  each  element  of  the  line  scan  camera  is  proportional  to  the  spectral  density 
Gd  (v)  of  the  combined  reference  and  sample  light,  which  can  be  expressed  as 

GAv)  =  Gs(v)\  +  ZRn  +  2Z  <jR„Rm  cos[2;rv(r„  - rm)]+  2E ^ cos[2^v(r„  -r,.)]j,  (l) 

v  n  n^m  n  3 

where  vis  the  light  frequency;  R„  is  the  normalized  intensity  reflection  representing  the  contribution  to  the 
collected  sample  light  by  the  nth  scatterer;  Gs  (v)  is  the  spectral  density  of  the  light  source;  the  reflection  of  the 

reference  arm  is  assumed  to  be  unity;  distances  are  represented  by  propagation  times  rn ,  rm  of  the  light 
reflected  by  the  nth  and  mth  scatterers  in  the  sample  and  zr  reflected  by  the  reference  mirror  in  the  reference 
arm  and  summation  is  across  all  axial  depths  in  the  sample  beam. 

The  spectral-domain  signal  can  be  transformed  to  the  time-domain  by  using  the  Wiener-Khinchin  theorem  [15]: 
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(2) 


r(r)  =<  u(t)u  *(t  +  T )  >=  FT-1  [G,  (v)] 

where  T(r)is  the  autocorrelation  function  of  the  source  light;  u(t  )  is  the  amplitude  of  the  electric  field  of  the 

light;  the  angle  brackets  denote  integration  over  time  and  FT~X  denotes  the  inverse  Fourier  transform.  By  taking 
an  inverse  Fourier  transformation  of  Eq.  (1),  we  obtain  the  time-domain  intensity  signal: 

7(r)  =  r(r)  +  r(r)Zi?„  +2E  ^Ar[r  +  2(r„  -  rm )]  +21  ^T[t  ±  2(z„ -zr)],  (3) 

n  n^m  n 

In  Eqs.  (1)  and  (3),  the  third  terms  are  the  mutual  interference  for  all  light  scattered  within  the  sample  expressed 
in  the  frequency  domain  and  the  time  domain,  respectively,  and  the  last  terms  contain  the  interference  between 
the  scattered  sample  light  and  the  reference  light  from  which  an  OCT  A-scan  is  calculated  [16]. 

The  discrete  Fourier  transformation  that  yields  Eq.  (3)  requires  even  sampling  in  v .  In  a  spectrometer,  however, 
the  spectrum  is  evenly  spread  along  the  wavelength  (here  we  do  not  consider  the  nonlinear  terms,  detailed 
analysis  of  the  wavelength  distribution  will  be  discussed  in  0)  and  the  acquired  raw  spectrum  must  be 
interpolated  to  get  the  correct  OCT  signal. 

Calculation  of  the  OCT fundus  image  [1 7] 

The  contrast  in  a  SLO  image  is  provided  by  the  lateral  distribution  of  X  R„  across  the  retina.  To  construct  a 
fundus  intensity  image  from  the  OCT  data  set  we  need  to  extract  the  intensity  term  X  Rn  that  is  contained  in 

both  the  second  and  fourth  terms  of  Eqs.  (1)  and  (3).  There  are  multiple  methods  for  extracting  the  intensity;  the 
method  selected  for  a  particular  application  will  depend  on  desired  speed  and  accuracy. 

One  method  is  to  use  the  non-interference  terms  in  the  frequency  domain  to  construct  the  intensity  image.  We 
noticed  that  the  cosine  terms  in  Eq.(l)  have  many  cycles  across  the  spectrum  and  will  sum  to  (approximately) 
zero,  leaving  only  the  constant  terms.  As  a  result,  when  we  sum  Eq.  (1)  across  v,  we  have 

Fvi(x,y)  =  G^1  +  X2?„j,  (4) 

where  Fv\ (x,  v)  is  the  output  of  the  processing  method  for  an  A-line  at  lateral  scan  point  (x,  y)  on  the  fundus  and 
Gs  is  the  total  source  power. 

Another  method  to  derive  the  fundus  intensity  is  to  use  the  fourth  term  of  Eq.  (1)  by  separating  the  oscillatory 
component  of  Eq.  (1)  from  the  slow  variation  and  recognizing  that,  for  retina  and  other  low  reflectance  samples, 
the  third  term  is  small  relative  to  the  fourth  term.  One  way  to  achieve  this  is  first  to  remove  the  low  frequency 
component  in  Eq.  (1)  by  high  pass  filtering  the  detected  spectrum,  then  squaring  the  remaining  oscillatory 
component  and  summing  over  the  spectrum.  The  result  can  be  expressed  as 

Fv 2  ( x ,  y)  =  X  kx  (v)  cos[2^v(rn  -  zr  )]}2  =  Ag)  X  Rn  ,  (5) 

where  Fv2  (x,  y)  is  the  intensity  calculated  in  the  frequency  domain,  which  can  be  displayed  directly  to  produce 
an  intensity  image.  According  to  Parseval’s  theorem,  Fv2(x,y)  =  Ft(x,y)  where  Fr(x,y)  is  the  intensity 
calculated  in  time  domain.  Ft  (x,  y)  can  be  acquired  from  the  calculated  OCT  signal  in  Eq.  (3)  by  squaring  and 
summing  the  values  at  all  axial  positions  except  those  near  z  =  0  .  We  have 

F, (x, y)  =  X  kx  ±  2 (rn  -  zr )]}2 .  (6) 

r  v  n  1 

Wavelength  distribution  in  the  spectrometer 

In  a  spectral-domain  OCT  system  the  accuracy  of  the  wavelength  distribution  in  the  spectrometer  can  severely 
affect  the  depth  resolution  of  the  system  and  the  signal  to  noise  ratio  (SNR).  Therefore,  calibration  of  the 
wavelength  distribution  in  the  spectrometer  is  one  important  step  for  achieving  depth  resolution  close  to  the 
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theoretical  prediction  [18].  Fig.  16  shows  the  simulation  results  on  the  effect  of  wavelength  distribution  error  in 
the  spectrometer  on  the  depth  resolution  in  spectral-domain  OCT.  Two  boundaries  with  different  reflectivity 
were  simulated  as  the  sample.  Fig.  16a  shows  the  situation  when  there  is  no  error  in  the  wavelength  distribution. 
Fig.  16b  shows  the  situation  when  the  nonlinear  terms  in  the  wavelength  distribution  exist  but  were  omitted. 
The  resulted  depth  dependent  resolution  distortion  can  be  clearly  seen  in  the  figure. 


Depth  (pixel) 


(a) 


(b) 


Fig.  16.  Simulation  of  the  effect  of  error  of  wavelength  distribution  in  the  spectrometer  on  the  depth 
resolution  in  spectral-domain  OCT.  Two  boundaries  with  different  reflectivity  were  simulated  as  the 
sample,  (a):  there  is  no  error  in  the  wavelength  distribution;  (b):  the  nonlinear  terms  in  the  wavelength 
distribution  exist  but  were  omitted.  We  can  see  clearly  the  depth  dependent  resolution  distortion  in  (b). 


The  theoretical  wavelength  distribution  in  the  spectrometer  can  be  worked  out  using  Fig.  17a.  The  grating 
equation  for  the  first  order  diffraction  can  be  expressed  as 


X  =  d  sin  6i  +  d  sin(#  +  #,.), 


(7) 


where  6i  is  the  incident  angle  for  the  transmission  grating;  6t+  6  is  the  diffraction  angle  for  wavelength  X;  d  is 
spacing  of  the  grooves  of  the  grating.  Assuming  that  the  center  wavelength  of  the  light  source  falls  on  the  center 
pixel  of  the  CCD  camera  the  wavelength  distribution  on  the  CCD  camera  can  be  expressed  as 


X  =  ^  + 
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=  X0  +  a1(xd -Vm  12 )-a2(xd -Vm  / 2)2 -a3(xd-Vm /2)3 
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where  Xo  is  the  center  wavelength  of  the  light  source;  x  is  the  distance  on  the  CCD  imaging  plane  from  the 
center; /is  focal  length  of  the  imaging  lens;  a\,  a2,  and  a?,  are  coefficients  for  the  first  order,  second  order,  and 
third  order  terms;  Xd  and  Vm  are  the  pixel  number  and  the  maximal  pixel  number  of  the  CCD  camera,  0<  xj  <  Vm. 
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Fig.  17.  (a)  Illustration  of  the  light  diffraction  in  the  spectrometer.  9t  is  the  incident  angle  of  the 
incoming  light.  The  configuration  of  the  system  meets  the  Littrow  condition.  0  is  the  diffraction  angle 
in  reference  to  the  center  wavelength,  (b)  An  example  for  the  spectromerter  calibration:  the  position 
for  each  atomic  line  predicted  by  the  theoretical  model  (*)  was  compared  with  the  measured  spectrum 
of  a  Neon  spectral  lamp. 

Fig.  17b  shows  one  example  of  the  spectrometer  calibration  using  a  Neon  spectral  calibration  lamp  (6032, 
Newport  Corporation),  where  the  position  of  each  atomic  line  predicted  by  the  theoretical  model  (*  in  the 
figure)  was  compared  with  the  measured  spectrum  of  the  lamp.  In  this  example  the  measured  position  of  each 
atomic  line  of  the  spectral  lamp  matched  the  theoretical  model  preferably.  When  errors  occur,  the  coefficients 
a\,  cii,  and  <23  need  to  be  tuned. 

To  further  eliminate  any  effect  from  the  possible  residual  error  in  the  wavelength  distribution  after  calibration 
with  the  spectral  lamp  a  method  reported  in  by  C.  Dorrer  [18]  will  be  employed.  Namely  we  can  use  a  mirror  in 
the  sample  arm  and  measure  the  spectral  interferogram  with  different  delays.  By  curve  fitting  on  the  phase 
differences  for  different  delay  differences  the  non-linear  term  in  the  wavelength  calibration  error  can  be 
extracted  and  then  eliminated  during  signal  processing  for  OCT  imaging. 

Dispersion  compensation 

Dispersion  is  caused  by  the  dependence  of  the  refractive  index  of  the  optical  medium  on  the  optical  frequency. 
In  a  dispersive  optical  medium  light  of  different  wavelengths  travels  at  different  speed.  Dispersion  mismatch 
between  the  sample  and  reference  arms  in  an  interferometer  causes  distortions  to  the  point  spread  function,  then 
the  depth  resolution,  of  the  system.  Dispersion  between  the  two  arms  in  an  interferometer  needs  to  be  carefully 
matched  in  the  entire  bandwidth  of  the  broadband  light  source  in  order  to  achieve  optimal  depth  resolution  in  an 
ultra-high  resolution  OCT  [19].  Dispersion  in  the  sample  arm  in  retinal  OCT  imaging  is  predominantly 
attributed  to  the  media  in  front  of  the  retina  including  the  ocular  tissue  and  optical  materials  in  the  sample  arm. 
As  a  result,  dispersion  variation  in  the  entire  imaging  range  is  negligible  [20]. 

In  spectral-domain  OCT  dispersion  mismatch  between  the  sample  and  reference  arms  can  be  compensated 
during  signal  processing.  The  phase  of  the  interferogram  for  each  axial  scan  can  be  expressed  as 

®(<z>)  =  O0  +  cx(a>-  coQ)  +  c2(a>-  a>0)2  +  c3(a>  -  a>0)3 ,  (9) 

where  a>0  is  the  center  angular  frequency  of  the  light  source.  The  first  order  term  (the  second  term)  is  the  OCT 
signal.  The  second  and  third  order  terms  (the  third  and  fourth  terms)  are  caused  by  the  second  and  third  order 
dispersion  mismatch  between  the  two  arms.  To  compensate  for  the  dispersion  mismatch  we  can  construct  a 
complex  compensation  phase  function: 

Gc (i co )  =  exp [-ic'2 (co  -co0)2  -  ic\ (co  -  co0 )3  ] .  (10) 

After  removing  the  DC  terms  in  the  detected  interferogram  expressed  in  Eq.  (1)  we  multiply  the  signal  with 
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(11) 


Gc  ( co )  and  we  have 

Gc(co)Gd(co)  =  exp[-ic'2(co-co0)2  -ic'3(co- co0)3]2'Zy[R^ cos[co(Tn  -  rr)  +  c2(a> -  a>0)2  +  c3(co  -  a>0) 

n 

=  X  exp[/ffl(TB  -Tr)  +  i(c2  -  c2 )(co -  co0 )2  +  (c3  -  c3 )(a> -  co0 )3 ] 

n 

+  X  exp[/<2?(r„  -  rr )  +  i(c2  +  c'2)(a> -  a>0)2  +  (c3  +  c3 )(®  -  co0 )3  ] 

n 


where  we  omitted  the  initial  phase.  When  c2  =  c2  andc3  =c3,  by  taking  the  inverse  Fourier  transformation  we 
have 


FT1  [Gc  (co)Gd  (®)]  =  FT 1  fe  ^  exp [m(rn  -  rr  )]j 

r  ^ 


+  FT  1  jl  y[R^ exp[/ft>(r„  -rr)  +  2 ic2 (co-a>0)2  +  2 ic3 {co -  a>0 ) 3  j 


(12) 


The  first  term  in  the  expression  is  the  dispersion  free  OCT  signal  while  the  second  term  is  the  mirror  image  in 
which  the  dispersion  is  doubled.  With  this  method  the  real  image  is  compensated  for  the  dispersion  mismatch 
while  the  mirror  image  is  blurred,  which  has  the  additional  advantage  for  discriminating  the  real  image  from  its 
mirror  for  the  operator. 

An  iterative  algorithm  similar  to  the  one  used  by  M.  Wojtkowski  et  al  [20]  for  finding  the  correct  coefficients 
for  c2  =c2  and  c3  =c3  will  be  used.  The  task  of  the  algorithm  is  to  find  the  two  coefficients  that  make  the 

energy  in  each  axial  scan  the  most  concentrated,  i.e.  the  image  is  the  sharpest.  The  flowchart  for  the  procedure 
of  dispersion  compensation  is  shown  in  Fig.  18.  In  the  procedure  the  last  zero  padding  is  to  ensure  that  after  the 
inverse  Fourier  transformation  the  pixel  spacing  is  smaller  than  half  of  the  designed  depth  resolution.  The 
algorithm  eliminated  the  process  for  calculating  the  phase  of  the  signal  thus  avoided  the  possible  unwrapping 
error. 


Fig.  18.  Flow  chart  for  the  iterative  dispersion  compensation  procedure.  The  iteration 
continues  until  the  OCT  signal  is  the  sharpest  measured  with  the  sharpness  metric.  Usually 
the  processes  for  finding  the  second  and  third  order  compensation  coefficients  are  separated. 
The  last  zero  padding  is  to  ensure  that  after  the  inverse  Fourier  transformation  the  pixel 
spacing  is  smaller  than  half  of  the  designed  depth  resolution. 


As  shown  in  Fig.  19,  the  calibrated  depth  resolution  is  3.8  pm  in  air  and  ~3pm  the  tissue,  which  was  corrected 
with  the  refractive  index  of  biological  tissues  (~1.4). 
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Fig.  19.  measured  spectrum  of  the  light  source  and  the  point-spread  function  of  the  ultra-high 
resolution  OCT. 

Test  the  OCT  system  on  normal  human  eye 

The  OCT  system  was  tested  on  imaging  normal  human  eye  (the  Pi’s  eye).  Fig.  20  shows  the  OCT  image  of  a 
normal  human  retina.  From  the  image  we  can  see  that  all  the  sub-retinal  layers  can  be  seen  clearly,  which 
demonstrated  the  resolving  capability  of  the  ultra-high  resolution  OCT  system. 


Fig.  20.  OCT  image  of  the  normal  human  retina. 

By  adding  switchable  optics  in  the  sample  arm,  the  OCT  system  was  capable  of  imaging  both  the  anterior 
segments  and  the  retina  of  the  eye.  Fig.  21  shows  the  acquired  images  of  the  cornea  and  conjunctiva  of  a  normal 
eye  with  contact  lenses.  By  adding  artificial  tears  and  after  several  blinks  the  OCT  system  can  not  only  reveal 
the  pre-lens  and  post-lens  tear  films  but  also  all  the  details  of  the  anatomy  of  the  cornea:  the  corneal  epithelium, 
basal  cell  layer,  Bowman’s  membrane,  stroma,  and  endothelium. 
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Fig.  21  OCT  images  of  the  anterior  segment  of  the  eye.  (a)  a  magnified  view  of  the  cornea 
with  soft  contact  lens;  (b)  image  of  the  conjunctiva  with  soft  contact  lens. 


KEY  RESEARCH  ACCOMPLISHMENTS 

■  Developed  methods  for  accurate  vertical  alignment  and  matching,  including  the  overlay  of  a  color  or 
red- free,  fundus  photograph  on  an  OCT  image  for  registration. 

■  Developed  the  mathematical  algorithms  to  analyze  and  quantify  OCT  datasets  automatically  including 
3D  segmentation  algorithm  for  quantitative  mapping  of  the  retinal  thickness  using  3D  OCT  data  for 
normal  and  diseased  human  eyes. 

■  Tested  the  accuracy  of  the  segmentation  algorithm  by  comparing  the  segmentation  results  done  by  the 
algorithm  and  by  experts,  respectively. 

■  Further  improved  the  resolution  of  the  spectral  OCT.  The  developed  OCT  system  has  a  depth  resolution 
of  ~3  pm  in  tissue  and  is  capable  of  resolving  all  the  sub-retinal  features  of  the  retina. 


REPORTABLE  OUTCOMES 
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calculation  using  spectral  domain  optical  coherence  tomography,"  Opt.  Express  15,  15193-15206  (2007). 
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1) .  S.  Jiao,  J.  Wang,  and  C.  A.  Puliafito,  “Contrast  Enhancement  for  Imaging  the  Tears  and  Contact  Lens  with 
Optical  Coherence  Tomography”,  ARVO  2008. 

2) .  H.M.  Wehbe,  M.  Ruggeri,  S.  Jiao,  C.  Puliafito,  “Quantitative  Retinal  Blood  Flow  Measurement  and 
Calibration  Using  Spectral  Domain  Optical  Coherence  Tomography”,  ARVO  2008. 

CONCLUSION 

High  resolution  OCT  images  can  reveal  the  detailed  anatomical  structures  of  both  the  anterior  and  posterior 
segments  of  the  eye.  Ultra-high  resolution  SD-OCT  provides  a  powerful  tool  for  visualizing  the  3D  microscopic 
structures  of  the  eye  in  vivo.  Incorporating  with  segmentation  algorithms  quantitative  information  like  the  3D 
geometry  of  the  retina,  retinal  thickness  map,  and  tear  volume  can  be  extracted  from  the  measured  OCT  images, 
which  make  possible  more  objective  quantitative  evaluation  of  eye  diseases.  It  will  also  simplify  and  accelerate 
diagnosis  and  monitoring  of  retinal  disease  by  enabling  prompt  triage  and  therapy  implementation  by 
optometrists  and  non-ophthalmologist  physicians.  The  development  of  registration  algorithms  makes  possible 
more  accurate  evaluation  of  images  taken  at  progressive  intervals,  which  will  better  enable  ophthalmologists  to 
follow  the  course  of  the  disease  and  its  treatment.  Furthermore,  improved  resolution  of  SD-OCT  will  enable 
telemedicine  consultation  on  retinal  disease  for  non-ophthalmologist  health  professionals,  potentially  preventing 
vision  loss  and  blindness.  Next  step  of  the  study  should  be  focused  on  optimization  of  the  system  hardware  and 
registration  software  to  make  it  more  stable  and  robust. 
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APPENDICES 


Automatic  retinal  blood  flow  calculation  using 
spectral  domain  optical  coherence  tomography 
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Abstract:  Optical  Doppler  tomography  (ODT)  is  a  branch  of  optical 
coherence  tomography  (OCT)  that  can  measure  the  speed  of  a  blood  flow 
by  measuring  the  Doppler  shift  impinged  on  the  probing  sample  light  by  the 
moving  blood  cells.  However,,  the  measured  speed  of  blood  flow  is  a 
function  of  the  Doppler  angle,  which  needs  to  be  determined  in  order  to 
calculate  the  absolute  velocity  of  the  blood  flow  inside  a  vessel.  We 
developed  a  technique  that  can  extract  the  Doppler  angle  from  the  3D  data 
measured  with  spectral-domain  OCT,  which  needs  to  extract  the  lateral  and 
depth  coordinates  of  a  vessel  in  each  measured  ODT  and  OCT  image.  The 
lateral  coordinates  and  the  diameter  of  a  blood  vessel  were  first  extracted  m 
each  OCT  structural  image  by  using  the  technique  of  blood  vessel 
shadowgram,  a  technique  first  developed  by  us  for  enhancing  the  retinal 
blood  vessel  contrast  in  the  ew  face  view  of  the  3D  OCT.  The  depth 
coordinate  of  a  vessel  was  then  determined  by  us  mg  a  circular  averaging 
filter  moving  m  the  depth  direction  along  the  axis  passing  through  the 
vessel  center  in  the  ODT  image.  The  Doppler  angle  was  then  calculated 
from  the  extracted  coordinates  of  the  blood  vessel.  The  technique  was 
applied  in  blood  flow  measurements  in  retinal  blood  vessels,,  which  has 
potential  impact  on  the  study  and  diagnosis  of  blinding  diseases  like 
glaucoma. 

©2007  Optica  l  Society  of  America 

QCIS  codes:  (110.4500)  Optical  coherence  tomography;  (120.3 890)  medical  optics 
instrumentation:  {170.45 B0)  optical  diagnostics  fo:  medicine. 
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1.  Introduction: 

Glaucoma  is  one  of  the  leading  causes  of  blindness  in  the  world  and  is  usually  associated  with 
increased  intraocular  pressure  (IOP).  There  are  two  theories  for  the  pathogenesis  of 
glaucomatous  optic  neuropathy  (GON):  the  mechanical  and  the  vascular  theory-;,  both  of  which 
have  been  the  debate  subject  of  multiple  research  groups  throughout  the  past  150  years.  [1,2- 
3]  The  vascular  theory  considers  GOX  to  be  the  result  of  ischemia  caused  by  the  elevated  IOP 
or  other  risk  factors  obstructing  the  blood  flow.  A  tool  providing  accurate  quantitative 
structural  and  blood  flow  information  will  benefit  the  study  of  the  etiology  of  glaucoma  as 
well  as  the  developments  of  new  therapie  s  by  monitoring  the  treatment  effects  non-invasively. 
Optical  Doppler  tomography”  (GDT)  is  a  branch  of  optical  coherence  tomography  (OCT)  that 
can  measure  the  spatially  resolved  speed  of  a  blood  flow  by  measuring  the  Doppler  shift 
caused  by  the  moving  blood  cells  to  the  probing  sample  light.  [4.5]  The  measured  Doppler 

shift  (/^  )  is  related  to  the  velocity  by: 


A  = 


cos# . 


(i) 


where  Va  is  the  absolute  velocity’  of  the  moving  blood  cells,  is  the  center  wavelength  of 

the  light  source,  n  is  the  refractive  index  of  the  sample  and  9  is  the  Doppler  angle.  From  the 
above  equation  it  is  clear  that  the  calculated  velocity  from  the  measured  Doppler  shift  is  the 
projection  of  the  absolute  velocity  on  the  direction  of  the  incident  probing  light.  Accordingly, 
in  order  to  calculate  the  absolute  velocity  from  the  measured  Doppler  shift  we  need  to  know 
the  Doppler  angle  9 .  However  9  is  not  only  an  unknown  parameter  in  retinal  ODT  imaging 
but  also  changes  at  different  imaging  time. 

Spectral-domain  optical  coherence  tomography  (SD-OCT)  is  a  recently  developed  high 
speed  OCT  technology  that  provides  high  resolution  three  dimensional  imaging  of  biological 
tissues.  SD-OCT  provides  the  means  of  calculating  the  Doppler  angle  from  the  acquired  3D 
data,  which  provides  structural  (3D  orientations  and  the  blood  vessel  diameters)  and  Doppler 
information  of  the  retinal  blood  vessels.  [6]  In  this  paper  we  report  on  our  study  on  the 
automatic  extraction  of  the  parameters  of  the  retinal  blood  vessels  from  the  images  acquired 
with  a  high  resolution  high  speed  SD-OCT.  The  retinal  blood  flow  can  be  calculated  upon 
acquisition  of  these  parameters  together  with  the  Doppler  information. 

2.  Materials  and  methods 

2.1  Ulfra  high  resolution  OCT  sy  stem 

A  high-speed  high  resolution  3D  SD-OCT  was  built  for  the  investigation.  In  the  SD-OCT 
system,  the  low- coherence  light  from  a  three-module  superluminescent  diode  (T-840 
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Broadlighter.  Superlum  Diodes  Ltd.  Moscow,  Russia)  with  center  wavelength  of  840  mu  and 
FWHM  bandwidth  of  100  nm  is  coupled  into  an  optical  fiber-based  Michelson  interferometer. 
In  the  sample  arm.  the  sample  light  is  delivered  to  the  retina  by  a  modified  optical  head  of  an 
OCT  2  system  (Carl  Zeiss  Meditec  Inc..  Dublin,  CA).  In  the  detection  arm.  a  spectrometer 
consisting  of  a  collimating  lens,  a  transmission  grating  (1200  lme/miu),  a  multi-element 
imaging  lens  (f  =  ISO  mm),  and  a  line  scan  CCD  camera  (Avhva-SM2-CL-2014.  2048  pixels 
with  14  micron  pixel  size  operating  in  12 -bit  mode)  was  used  to  detect  the  combined 
reference  and  sample  light.  The  calculated  spectral  resolution  of  the  spectrometer  was  0.055 
mil.  which  corresponds  to  a  detectable  imaging  depth  range  of  3.1  mm  in  air.  The  calibrated 
depth  resolution  in  the  tissue  is  -3|im.  which  was  corrected  with  the  refractive  index  of 
biological  tissues.  The  power  of  the  sample  light  was  lowered  to  7.5 OllW  by  adjusting  the 
source  power  to  ensure  that  the  light  intensity  delivered  to  the  eye  was  within  the  ANSI 
standard.  OCT  scans  consisting  of  a  total  number  of  65536  depth  scans  (A-scans)  takes  2.7 
seconds. 

2.2  Doppler  imaging 

The  measured  raw  spectral  interference  signals  were  processed  by  using  the  standard 
algorithm  for  spectral  domain  OCT  to  get  the  complex  signal  [  V(t)  ]  in  the  tune  domain.  To 
calculate  the  Doppler  shift  we  first  calculated  the  product  of  the  complex  signals  for  the 
adjacent  A- lines: 


r1+1(0  *  r*  (o  -  w-*(,a  (2) 

where  I”*  I/)  is  the  complex  conjugate  of  T(r) :  A  is  the  amplitude  of  the  product;  (pi  is  the 
phase  of  l~(r);  r  is  the  A-line  number.  The  phase  difference  among  the  adjacent  A- 
lines  A(Pj  =  (pi+l  -<Pj  can  then  be  calculated.  Calculating  the  phase  difference  with  Eq.  (2)  has 
the  advantage  of  avoiding  the  problem  of  phase  unwrapping. 

The  projected  flow  speed  on  the  direction  of  the  incident  sample  light  can  be  calculated 
from  A  (Pj : 


—K  <  Apj  <x  , 


_  A-line 

4m? 

_Mi  -line  ^  hf.i  -j'rpig 

4n 


(3) 


where  v  is  the  projection  of  the  absolute  velocity  Vff  along  the  depth  direction,  Aq  is  the 

center  wavelength  of  the  light  source  (840  nm).  ne  is  the  axial  scan  frequency  (24000  A- 
lines/sec),  n  is  the  index  of  refraction  of  the  sample  (—14)  Accordingly,  the  maximal  speed 
that  can  be  detected  by  our  system  without  phase  wrapping  is  ±3.6  mm/sec. 

The  Doppler  image  contains  bulk  motion  artifacts  that  is  additive  to  the  Doppler  shift 
induced  by  the  blood  flow.  Bulk  motion  effect  can  be  quantified  for  each  A-line  by  using  the 
histogram  technique  [7].  Briefly  speaking,  the  histogram  H[Atp})  of  the  phase  difference  A 
was  calculated  along  each  individual  line  of  the  Doppler  image.  H{A(p-  ]  contained  N  bins 

varying  in  the  interval  [-/r,jr]  with  a  step  size  of  A ^  f  n:tm  ?  where  and 

Affy,-r,  are  the  maximum  and  minimum  of  A$r ,  respectively.  The  phase  shift  A<pbulk  caused 
by  bulk  motion  was  obtained  by  locating  the  peak  of  the  histogram 

=  max[tf(A^ )],  -n  <  H{ )<k  (4) 
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A<pbuIk  is  subtracted  from  its  corresponding  A(p.  to  compensate  for  the  bulk  motion  artifact 
for  each  line  of  the  ODT  image. 

2.3  Scan  patterns  and  eve  movement  compensation 

The  first  scan  pattern  we  used  for  the  study  is  a  set  of  evenly  spaced  concentric  circles  around 
the  optic  disc  of  the  eye.  This  scan  pattern  was  chosen  to  cover  all  major  retinal  arteries  and 
veins.  This  scan  pattern  was  used  to  test  the  algorithms  for  blood  vessel  detection.  As 
illustrated  in  Fig.  1(a),  the  circular  scan  pattern  consisted  of  56  normal  density  (1024)  and  one 
high  density  (8192)  depth  scans  (A-lines).  The  radius  of  the  circles  spanned  evenly  from  0.73 
to  2.73mm  on  the  fundus.  The  high  density  scan  had  a  scan  radius  of  1.73mm.  which  is  the 
same  as  that  used  in  the  commercial  time-domain  OCT  machine  for  glaucoma  imaging. [8 ]r  A 
donut  shaped  en  face  OCT  fundus  image  [11]  can  be  generated  and  registered  on  the 
corresponding  fundus  photograph  for  the  same  eye.  By  comparing  the  OCT  fundus  image  and 
the  fundus  photograph  arteries  and  veins  can  be  recognized.  An  ODT  image  can  be  generated 
from  the  high  density  scan.  Hie  ODT  image  helps  not  only  verify  the  accuracy  of  the 
calculation  of  the  lateral  coordinates  of  each  recognizable  blood  vessel  but  also  locate  the 
depth  coordinate  of  each  blood  vessel  on  the  high  density  OCT  unage. 


Fig.  1.  (a).  Circular  sea  a  pattern  centered  on  tbe  optLC  disk.  Tbe  black 
circles  represent  die  norma]  den&Lty  scans  (1024  A-lines):  die  white  arcle 
represent  the  high  density  scan{S192  A-lines).  (b)Arc  ^ can  pattern.  Scans 
1,2-3.62.63  are  used  to  calculate  die  DoppLer  angle.  Scan  64  is  used  for 
alignment.  Scans  4-61  were  scanning  the  same  area. 

Another  scan  pattern  we  used  is  arc  shaped  [Fig.  1(b)]  for  imaging  individual  blood 
vessels.  The  scan  pattern  consisted  of  63  concentric  arcs,  each  of  which  has  1024  A-lines, 
subtending  an  angle  of  M  /4  with  a  radius  spanning  from  1.48  to  1.98mm.  and  one  depth 
alignment  scan  (reference  scan,  scan  No.  64)  for  eye  movement  compensation.  58  scans  (scan 
No.4  to  No. 61)  were  repeated  at  a  radius  of  1.73  mm.  Tins  arrangement  of  the  scans  is  used  to 
calculate  the  Doppler  angle  and  the  flow  dynamics.  A  Doppler  image  was  calculated  for  each 
of  the  arc  scans,  therefore  allowing  accurate  3D  coordinate  calculations. 

Also  illustrated  in  Fig.  1(a)  is  the  coordinate  system  we  used  for  the  entire  study,  where 
the  X-axis  represents  the  horizontal  position,  the  Y-axis  represents  the  vertical  position,  and 
the  Z-axis  represents  the  depth  position.  We  also  defined  r  (1.48mm  <  r  <  l_98mm  )  as  the 
scan  radius  and  a  (0  <a<  2/r)  as  the  scan  angle,  both  can  be  pre -determined  by  the  scan 
data.  For  each  pixel  on  the  circular  OCT  image  the  lateral  coordinates  (x,  y)  can  be  calculated 
as  x  =  r  cos  a  and  y  =  r  sm  a . 
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Involuntary  eye  movements  during  image  acquisition  can  cause  distortions  to  the  retinal 
OCT  images  a?  well  as  the  location?  of  retinal  vessels  due  to  limited  imaging  speed,  which  in 
turn  add  an  error  to  the  Doppler  angle  calculation.  The  ultimate  solution  to  this  problem  is  to 
increase  the  imaging  speed.  With  current  commercially  available  camera?  for  the 
spectrometer,  the  space  for  improving  imaging  speed  is  very  limited.  As  a  result, 
compensation  for  the  distortion?  caused  by  eye  movements  with  post-processing  is  necessary. 
To  compensate  for  the  eye  movements  a  linear  fast  reference  scan  crossing  all  the  arc  scans 
was  added  at  the  end  of  the  arc  scan.  During  the  fast  linear  scan  eye  movement  can  usually  be 
neglected,  the  linear  scan  provides  a  reference  for  all  the  arc  scan  images  in  the  depth 
direction.  The  depth  position  for  each  arc  scan  can  be  adjusted  according  to  this  reference 
scan  and  thus  provides  compensation  in  the  depth  direction.  Compensation  according  to  the 
reference  scan  unage  can  be  done  in  different  ways.  In  one  method  we  can  detect  the  front 
surface  of  all  the  images  and  each  B-scan  will  be  shifted  in  the  Z  direction  according  to  the 
difference  of  the  z  coordinates  of  the  surfaces  between  the  B-scan  and  the  reference  image  at 
the  corresponding  (x:  y)  position.  We  can  also  calculate  the  shift  for  each  B-scan  by  means  of 
the  correlation  coefficients  between  the  corresponding  A-lines  in  the  circular  and  reference  B- 
scans. 

The  blood  vessels  can  be  treated  as  straight  lines  in  the  small  arc  scan  area.  The  Doppler 
angle  can  be  calculated  using  the  following  relationship: 


cos$  = 


Ax2  +Av3  +  Az2 


(5) 


Where  Ax, ,  Ay  T  and  Ac  are  the  projections  of  the  vessel  segment  in  the  scanned  area  on  the  Xr 
Y,  and  Z  axes.  By  differentiation  of  Eq.  (5)  we  have 


d(cos&)  = 


-  Ax  Aid  (Ax)  -  Ay  Aid  (Ay)  +  (Av"  +  Ay")d(Ai) 
(Ax2  +  Ay2  +  A z2)3'2 


(6) 


In  retina  imaging,  #is  close  to  90°  and  Az  «  Ay,  Ay  .  As  a  result,  we  can  see  from  Eq.  (6)  the 

variation  of  the  Doppler  angle  is  more  sensitive  to  the  variation  of  Az  .  So  that  compensation 
for  eye  movement  in  the  z  direction  is  the  most  important  m  calculatmg  the  Doppler  angle. 


Fig  2.  (a),  virtual  B-scan  extracted  from  the  3D  data  at  the  location  of 
the  reference  scan.  The  red  line  shows  the  segmented  ELM  {inner 
limiting  membrane)  that  is  used  for  augmnen:.  (b)  The  reference  scan 
hmaae. 

Fig.  2  shows  the  extracted  image  (virtual  B-scan)  from  the  acquired  3D  data 
corresponding  to  the  positions  of  the  linear  reference  scan  together  with  the  reference  scan 
miage.  The  inner  limiting  membrane  (ELM)  can  be  extracted  by  using  our  segmentation 
algorithm  for  both  unages.  By  comparing  the  Z  coordinates  at  the  same  locations  of  the 
segmented  ELM  in  both  images,  a  compensation  curve  can  be  generated. 
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2.4  Automatic  quantitative  blood  vessel  detection 


2,4.1  Automatic  detection  of  lateral  coordinates 


To  calculate  the  Doppler  angle,  the  coordinates  of  the  blood  vessel  need  to  be  determined  on 
each  OCT  cross  sectional  image.  Our  strategy  is  to  first  determine  the  lateral  position  and  the 
diameter  of  each  blood  vessel  on  the  cross  sectional  OCT  intensity  image  and  then  determine 
the  depth  location  of  the  vessel  on  the  ODT  image.  In  an  OCT  cross  sectional  image  a  blood 
vessel  casts  a  shadow  behind  it,  which  exhibits  as  a  low  reflection  region  in  the  retinal  fundus 
reflection  graph.  However,  surface  reflections  can  cause  distortions  to  the  blood  vessel  profile 
and  add  difficulty  in  blood  vessel  automatic  detection.  Other  techniques  were  introduced  to 
minimize  the  effect  of  surface  reflection  in  retinal  blood  vessel  detection.  [12]  Here  we  used 
our  previously  published  method  known  as  the  shadowgram  [11]  to  improve  the  retinal  blood 
vessel  profile  and  minimize  the  surface  reflection  effect.  Using  this  technique  we  can  generate 
a  high  contrast  retinal  reflection  distribution  called  shadowgraph,  where  each  vessel  was 
characterized  as  a  low  reflection  region. 


Fig.  3.  Improved  blood  vessel  profile  urias  blood  vessel  -shadowgraph. 

(a)  Orisiaal  B-scan  image,  (b)  The  B-scan  image  after  the  -surface  layer 
was  removed.  The  calculaoed  fundus  reflection  distribution 
corresponding  to  the  B-^can  a  ad  tbe  fundus  shadowgraph  were 
superimposed  on  the  images. 

The  shadowgraph  was  generated  by  first  removing  the  surface  layers  of  the  retinal  OCT 
image  and  then  summing  all  the  pixel  intensity  values  along  each  A- sc  an.  Fig.  3  shows  the 
original  OCT  image,  the  image  after  the  surface  layer  was  removed,  and  the  calculated 
reflection  distributions  (the  vessel  profiles).  From  the  calculated  fundus  reflections  we  can  see 
that  not  only  the  vessel  contrast  but  also  the  vessel  profile  was  unproved  by  using  the 
shadowgraph,  winch  promises  a  more  accurate  detection  of  the  vessel  location  and  its 
diameter. 

To  extract  the  lateral  coordinates  the  shadowgraph  was  first  smoothed  by  using  a  smooth 
filter  (Savitzky-Golay,  length=3s.  weighing  factor=2.1  in  Matlab).  The  background  was 
removed  from  the  smoothed  curve  by  subtracting  the  low  pass  filtered  data.  After 
normalization,  thresholding  according  to  a  predetermined  value  was  applied.  The  blood  vessel 
locations  and  diameters  were  then  determined. 

2.4.2  Automatic  detection  of  Depth  coordinates 

Automatically  detecting  the  depth  position  of  a  blood  vessel  is  much  more  challenging  than 
detecting  the  lateral  coordinates.  Various  methods  may  be  used  on  either  the  intensity  image 
or  the  ODT  image.  We  first  reported  on  the  successful  detection  of  the  depth  position  of  a 
blood  vessel  on  an  ODT  image  after  the  lateral  coordinates  and  vessel  diameter  were 
determined.  [6]  The  same  method  was  used  in  the  current  study. 

To  detect  the  z  coordinate  of  the  center  of  a  blood  vessel,  comprehensive  information 
about  the  blood  vessel  was  used,,  i.e.  the  lateral  coordinates,  the  vessel  diameter,  and  the 
Doppler  shift  caused  by  the  blood  flow.  The  basic  idea  is  that  outside  the  blood  vessel  the 
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calculated  Doppler  shift  consists  of  only  random  noise  where  inside  the  vessel  the  Doppler 
shifts  are  correlated.  Accordingly,  a  circular  window  filter  moving  in  the  z  direction  along  the 
detected  vessel  position  can  be  constructed.  The  diameter  of  the  circular  window  should  be 
equal  or  less  than  the  diameter  of  the  blood  vessel.  The  filter  can  be  either  type  of  a  function 
that  examines  the  correlation  of  the  values  inside  the  widow.  The  simplest  type  of  the  function 
is  averaging. 

We  used  an  averaging  filter  for  the  detection  of  the  z  coordinate  of  a  blood  vessel.  An 
average  versus  z  curve  was  obtained  for  each  blood  vessel.  A  maximum  or  minimum  will  be 
reached,  depending  on  the  direction  of  the  blood  flow,  when  the  center  of  the  window 
coincides  with  the  center  of  the  blood  vessel  Therefore,  by  locating  the  maximum  or 
minimum  of  the  average  versus  z  curve  the  z  coordinates  of  the  center  of  a  blood  vessel  can 
be  determined. 

Fig.  4  shows  an  example  of  the  ODT  image  of  a  retinal  blood  vessel  in  an  arc  scan  around 
the  optic  disc  (a),  the  center  of  the  blood  vessel  marked  according  to  visual  measurement  (b), 
and  the  result  of  the  moving  circular  window  filtering.  In  this  case,  the  location  of  the 
minimum  corresponds  to  the  center  of  the  blood  vessel.  We  can  see  from  the  figure  that  the 
moving  circular  window  filtering  worked  well  in  locating  the  depth  position  of  the  vessel 
center. 


Fig.  4.  Depth  coordinate  detection  by  unng  a  circular  window  averaging 
(a)  The  original  ODT  image  of  an  arc  scan;  (b)  The  blood  vessel  center 
was  marked  by  VLmal  measurement;  (c)  Averaging  curve  along  the 
direction  of  the  A-line  passing  through  the  center  of  the  blood  vessel  with 
a  diameter  of -50  pixels.  Notice  die  position  of  the  peak  corresponds  with 
the  depth  position  of  the  center  of  the  blood  vessel. 


Upon  determination  of  the  coordinates  of  a  blood  vessel  the  Doppler  angle  can  be  calculated 
by  using  Eq.  (5).  The  absolute  velocity  of  the  blood  flow  inside  the  vessel  can  then  be 
calculated.  Knowing  the  absolute  average  velocity  (vff,  averaged  across  the  vessel)  of  the 
blood  flow  and  the  vessel  diameter  (d),  the  blood  flow  rate  (R)  can  be  calculated  as 

R=va/d2/  4.  (7) 
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Fig.  5  illustrates  the  procedure  for  the  data  processing 


Fig.  5.  Illustration  of  the  procedure  for  calculation  of  die  Doppler  angle  and  the  blood  flow. 


3.  Results  and  discussion 


Fig.  6.  Color  fundus  photograph  of  a  normal  human  eye  and  die 
corresponding  OCI  fundus  image  generated  from  the  circular  scans 
around  the  optic  disc. 

The  algorithms  for  detecting  the  vessel  coordinates  and  diameters  were  applied  to  the 
measured  data.  Fig.  6  shows  the  reconstructed  OCT  fundus  image  from  the  circular  scan 
around  the  optic  disc  and  the  corresponding  color  fundus  photograph  of  a  normal  eye  of  a 
volunteer.  By  comparing  both  images  retinal  arteries  and  veins  can  be  recognized  on  the  OCT 
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fundus  image.  The  OCT  fundus  image  was  used  to  determine  the  accuracy  of  the  algorithm 
for  the  detection  of  the  lateral  position  of  the  blood  vessels.  The  OCT  fundus  image  registered 
well  with  the  fundus  photograph.  Fig.  7  shows  the  step-by-step  data  processing  results  and  the 
detected  vessel  centers  and  vessel  boundaries  on  the  circular  OCT  and  the  ODI  images. 


Fl2.  7.  (a).  OCI  image  of  a  high  densicy  circular  scan  {8192  A-line^) 
around  tbe  opdc  disc  of  a  normal  human  eye:  (b).  The  OCT  image  after 
removal  of  the  surface  layers:  (c).  Tbe  original  shadowgraph:  (d).  The 
shadowgraph  after  background  correction  and  norma Sizarion;  (e).  tbe 
shadowgraph  after  thresholding:  (f).  Recognized  blood  vessel  centers 
and  boundaries  are  marked  on  the  OCT  Linage:  (g).  The  ODT  image  for 
the  same  OCT  scan:  (h).  Magnified  view  of  the  region  marked  in  (g). 
where  the  calculated  blood  vessel  centers  are  marked. 

To  test  the  accuracy  of  the  algorithm  the  eyes  of  four  normal  volunteers  were  imaged  and 
analyzed.  We  define  the  accuracy  for  the  detection  of  the  lateral  coordinates  as  the  percentage 
of  the  number  of  blood  vessels  automatically  detected  to  the  number  of  blood  vessels  detected 
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visually  on  the  circular  OCT  fundus  linage.  We  achieved  100%  accuracy  for  the  eyes  we  have 
imaged.  The  accuracy  for  the  detection  of  the  depth  coordinate  was  defined  by  comparing  the 
visually  determined  blood  vessel  center  on  the  Doppler  image  with  the  center  determined  by 
the  algorithm  (see  Fig.  8  as  an  example).  If  the  distance  between  the  two  centers  was  less  than 
half  of  the  radius  of  the  vessel,  we  define  it  as  a  success.  For  the  imaged  normal  eyes  we 
achieved  an  accuracy  of  84.4%  for  all  the  vessels.  For  large  blood  vessels  (diameter  larger 
than  50  pixels)  depth  positions  were  detected  with  an  accuracy  of  93.3%. 


Fig.  S.  (a),  the  detected  vessel  location  and  boundaries  in  an  are  scan 
image:  (b)  comparing  the  automatically  detected  vessel  center  (solid 
circle)  and  die  visually  detected  center  (open  circle)  on  the  ODT  image 
of  a  vessel. 

The  detected  blood  vessel  coordinates  contain  the  effect  of  eye  movement.  Before  the 
vessel  coordinates  were  used  for  calculation  of  the  Doppler  angle  eye  movement 
compensation  was  applied  to  the  arc  scans.  After  eye  movement  compensation  linear  fitting 
was  used  for  the  x,  y,  and  z  coordinates  of  the  blood  vessels.  The  angle  of  the  fitted  line  in 
reference  to  the  Z  axis  is  the  calculated  Doppler  angle.  Fig.  9  shows  the  process  and  results 
for  eye  movement  compensation  for  a  blood  vessel. 

After  calculation  of  the  Doppler  angle,  the  absolute  blood  flow  velocity  can  be  calculated. 
After  the  ODT  image  was  median  filtered  the  flow  speed  inside  the  blood  vessel  was  averaged 
across  the  vessel  for  each  tune  point  to  get  the  mean  of  the  flow  speed.  One  artery7  marked  on 
Fig.  10(a)  for  a  normal  eye  was  studied.  The  calculated  Doppler  angle  for  this  vessel  is  87.8°. 
The  calculated  absolute  velocity7  of  the  blood  flow  over  time  is  shown  in  the  curve  of  Fig. 
10(b).  The  blood  flow  velocity  over  the  time  of  measurement  was  calculated  to  be  30.4=9.5 
mm/s  (mean  and  standard  deviation).  The  standard  deviation  reflects  pulsating  behavior  of  the 
artery7.  The  volunteer’s  pulse  rate  measured  separately  immediately  after  the  OCT  imaging 
was  58  pulses/mm,  which  agrees  with  the  velocity7  calculation  in  Fig.  10b.  The  diameter  of  the 
vessel  in  the  measured  region  is  109  pm.  As  a  result,  the  average  blood  flow  in  this  vessel  is 
17  pl/min. 

One  vein  marked  on  the  fundus  photograph  shown  in  Fig.  11(a)  for  another  normal  eye 
was  also  studied.  The  calculated  Doppler  angle  for  this  vessel  is  87D.  The  calculated  absolute 
flow  velocity'  of  the  blood  flow,  averaged  across  the  vessel,  is  shown  in  the  in  Fig.  1 1(b).  The 
averaged  blood  flow  velocity  over  the  time  of  measurement  was  calculated  to  be  16.4=3.9 
mm/s  (mean  and  standard  deviation).  The  diameter  of  the  vessel  in  the  measured  region  is  51 
pm.  As  a  result,  the  average  blood  flow  in  this  vessel  is  2.01  pl/min. 
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00  (b) 


Fig.  '9.  The  detected  x.  y.  and  z  coordinate*  of  the  center  of  a  blood 
vessel  versus  the  scan  radms  r.  Linear  fitting  were  used  for 
compensating  variations  of  the  x  and  y  coordinates.  Tbe  coordinates  at 
i=1.73  nun  are  die  average  of  die  results  of  the  5£  repealed  arc  scan^. 

(a)  x  versus  r:  (b)  y  versus  r;  (c)  z  versus  r  before  and  after  eye 
movement  compensation.  Also  shown  in  (c)  is  the  range  of  the  z 
coordinates  of  tbe  vessel  center  for  the  53  repeated  arc  scans. 

ODT  images  for  the  artery  and  vein  are  also  shown  in  Fig.  10  and  Fig.  11.  During  the 
alignment  for  imaging  acquisition  a  real-time  ODT  image  was  displayed  for  the  position  at 
the  repeated  arc  scan.  The  ODT  image  for  the  blood  vessel  of  interest  was  optimized  by 
adjusting  the  optical  head  of  the  OCT  system,  which  is  equivalent  to  adjusting  the  Doppler 
angle  for  the  specific  vessel.  When  the  ODT  image  for  one  vessel  was  optimized  the  quality 
of  the  images  for  the  other  vessels  also  covered  by  the  scan  maybe  deteriorated.  We  speculate 
the  reason  for  this  deterioration  is  caused  by  that  the  Doppler  angle  was  closer  to  3GD. 
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Fig.  10  Tlie  test  result  of  an  artery  for  a  aona^  human  eye.  (a)  The 
fuadiv  photograph,  the  cross  ^ectLcnal  OCT  image  at  the  position 
marked  on  the  fucidiu  photograph,  aad  the  ODT  image,  (b)  The 
calculated  absolute  flow  velocity  averaged  across  the  vessel  area. 

The  calculated  absolute  velocity  and  flow  rates  for  the  artery  and  vein  compared  well  with 
results  obtained  using  other  technologies!  13]  The  tests  of  the  technique  on  current  limited 
number  of  eyes  were  successful  although  more  experiments  are  needed  to  test  the  accuracy 
and  repeatability.  Because  the  Doppler  angles  of  retinal  blood  vessels  are  close  to  90°.  the 
error  of  the  calculated  absolute  blood  flow  velocity  is  very  sensitive  to  the  error  of  the 
calculated  Doppler  angle.  As  analyzed  in  section  2.3  the  Doppler  angle  is  more  sensitive  to 
the  eye  movement  in  Z  direction  than  that  in  X  and  Y  directions,  which  is  important 
considering  that  there  is  also  no  effective  technique  in  the  compensation  for  the  movements  in 
the  X  and  Y  directions.  As  a  result,  in  our  current  technique  eye  movements  in  the  X  and  Y 
directions  were  corrected  only  with  linear  fitting. 

To  test  the  accuracy  of  the  calculation  of  the  Doppler  angle  and  the  blood  flow  rate  we 
measured  the  total  flow  rate  of  a  vein  before  and  after  a  bifurcation.  Two  measurements  at  the 
locations  shown  in  Fig.  12  were  taken.  For  the  first  measurement  we  scanned  the  area  that 
contained  the  vein  before  bifurcation  (parent  vessel),  while  in  the  second  measurement  we 
scanned  the  area  that  contained  both  branches  after  bifurcation  (daughter  vessels).  The  time 
interval  between  the  two  measurements  was  about  15  minutes,  during  which  the  subject  kept 
the  same  body  position.  The  two  measurements  were  taken  at  the  same  conditions  where  the 
room  was  kept  dark  to  minimize  the  outside  influence  to  the  blood  flow.  The  results  of  the 
measurements  are  shown  in  Table  1.  The  calculated  average  blood  flow  entering  the 
bifurcation  from  the  two  branches  was  2.23  plmin  while  the  average  blood  flow  leaving  the 
bifurcation  was  2.17  uTrnin.  The  result  provided  good  validation  of  our  technique  for  the 
calculation  of  the  Doppler  angle  and  the  blood  flow  rate.  We  are  planning  to  build  a  phantom 
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simulating  a  blood  vessel  with  adjustable  3D  orientation  and  flow  rate  to  further  test  the 
accuracy  of  our  algorithm. 


(b)  time  (sec) 

Fig.  11.  The  test  result  of  a  vein  for  a  normal  human  eye.  (a)  The  fundus 
photograph,  the  OCT  cross  sectional  image  at  the  position  marked  on  the 
fundus  photograph:  and  the  ODT  image.  <»  The  calculated  absolute  flow 
velocity  averaged  across  tbe  vessel  area. 


Fig.  12.  Color  fundus  photograph  of  a  normal  human  eye  with  markers 
indicating  the  location  of  the  scan  areas.  Vessel  1  corresponds  to  the  vein 
before  bifurcation,  while  vessels  2  and  3  represent  :he  vessel  branches  afeer 
bifurcation. 
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Table  1.  Calculated  parameters  of  die  vessels  shown  in  Fig.  12. 


Vessel  =? 

l 

2 

3 

Duineta:  (|on) 

103 

66 

54 

Velocity  (mm/ sec) 

4.3 

8 

4.12 

Doppler  Angle 

81.6* 

76.3d 

Flow  (fd-'min) 

2.17±0.57 

1.65+0.23 

0.58=0.07 

4.  Conclusion 

Doppler  angle  of  retinal  blood  vessels  including  arteries  and  veins  were  successfully 
calculated  by  using  the  comprehensive  information  provided  by  high  speed  SD-OCT — 
structural  information  from  the  OCT  intensity  image  and  speed  information  from  the  ODT 
image.  The  lateral  coordinates  of  a  blood  vessel  can  be  extracted  accurately  by  using  the 
technique  of  blood  vessel  shadowgram,  which  not  only  enhanced  contrast  of  the  blood  vessel 
against  the  reflecting  background  but  also  unproved  the  vessels  profile.  The  depth  coordinate 
of  a  blood  vessel  was  calculated  by  using  a  moving  circular  window  filter  in  the  ODT  image 
after  the  lateral  coordinates  and  the  vessel  diameter  were  extracted.  By  calculating  the 
Doppler  angle  of  a  blood  vessel  the  absolute  blood  flow  velocity'  and  the  blood  flow  rate  and 
be  calculated.  The  technique  was  successfully  tested  on  retinal  arteries  and  veins  for  normal 
human  eyes. 
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